﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace PIF1006_TP1
{

    public class EquationSystem
    {
        #region Attributes

        private Matrix matrixA;
        private Matrix matrixB;
        private Matrix matrixXResolved;

        #endregion

        public static EquationSystem Create(Matrix matrixA, Matrix matrixB)
        {
            EquationSystem equationSystem;

            try
            {
                if (matrixA.IsSquare)
                {
                    if ((matrixB.SizeY == 1 && matrixB.SizeX == matrixA.SizeX) ||
                        (matrixB.SizeX == 1 && matrixB.SizeY == matrixA.SizeY))
                    {
                        //Si la matrice B n'est pas dans le bon sens, on la transpose
                        if (matrixB.SizeX > matrixB.SizeY)
                            matrixB = matrixB.TransposedMatrix;
                        equationSystem = new EquationSystem();
                        equationSystem.matrixA = matrixA;
                        equationSystem.matrixB = matrixB;
                    }
                    else
                        throw new NotSameDimensionException("La taille de la matrice B n'est pas compatible avec la taille de la matrice A.");
                }
                else
                    throw new NotSquareMatrixException("La matrice A du système d'équation doit être carrée.");
            }
            catch (Exception exception)
            {
                throw exception;
            }
            return equationSystem;
        }

        private EquationSystem()
        {

        }

        #region Propreties

        public Matrix MatrixA
        {
            get { return matrixA; }
        }

        public Matrix MatrixB
        {
            get { return matrixB; }
        }

        public Matrix MatrixX
        {
            get { return matrixXResolved; }
        }

        public bool IsResolved
        {
            get { return (matrixXResolved != null); }
        }

        #endregion

        #region Methods

        public Matrix FindXUsingCramer()
        {
            Matrix matrixAi;
            double determinantA = matrixA.Determinant();

            if (determinantA != 0)
            {
                matrixAi = new Matrix(matrixB.SizeY, matrixB.SizeX);
                matrixXResolved = new Matrix(matrixA.SizeX, 1);

                //Pour chacune des colonnes de la matrice A
                for (int ctrCol = 0; ctrCol < matrixA.SizeX; ctrCol++)
                {
                    //Création d'une copie de la matrice A pour ne pas la modifier
                    matrixAi = matrixA.DeepCopy();
                    //On remplace une colonne par la colonne de la matrice B
                    for (int ctrLine = 0; ctrLine < matrixAi.SizeY; ctrLine++)
                        matrixAi[ctrLine, ctrCol] = matrixB[ctrLine, 0];

                    matrixXResolved[ctrCol, 0] = matrixAi.Determinant() / determinantA;
                }
            }
            else
                throw new Exception("Le déterminant doit être différent de 0");

            return matrixXResolved;
        }

        public Matrix FindXUsingMatrixInversion()
        {
            double determinantA = matrixA.Determinant();
            if (determinantA != 0)
                matrixXResolved = (matrixA.AdjugateMatrix.TransposedMatrix * (1 / determinantA)) * matrixB;
            else
                throw new Exception("Le déterminant doit être différent de 0");
            return matrixXResolved;
        }

        public Matrix FindXUsingJacobi(double epsilon, bool divergent = false, Matrix matrixX = null, Matrix matrixXDifference = null)
        {
            double largestDifferenceEpsilon = 0;
            double newDifferenceEpsilon = 0;
            Matrix subMatrixA;
            Matrix subMatrixX;
            Matrix matrixXTemp;
            double sumElements;
            double calcul;

            //Instanciation de la matrice de résultat au premier incrément
            if (matrixX == null)
                matrixX = new Matrix(matrixB.SizeY, matrixB.SizeX);
            matrixXTemp = matrixX.DeepCopy();
            //Instanciation du vecteur d'écart au premier incrément
            if (matrixXDifference == null)
                matrixXDifference = new Matrix(matrixB.SizeY, matrixB.SizeX);

            for (int ctrLineB = 0; ctrLineB < matrixB.SizeY; ctrLineB++)
            {
                subMatrixA = matrixA.SubMatrix(-1, ctrLineB);
                subMatrixX = matrixX.SubMatrix(ctrLineB, -1);
                sumElements = 0;
                calcul = 0;
                //On fait la somme des sous-matrices en éliminant la ligne et la colonne de l'élément courant
                for (int ctrColA = 0; ctrColA < subMatrixA.SizeX; ctrColA++)
                {
                    sumElements += subMatrixA[ctrLineB, ctrColA] * subMatrixX[ctrColA, 0];
                }
                calcul = (matrixB[ctrLineB, 0] - sumElements) / matrixA[ctrLineB, ctrLineB];

                //On calcule la nouvelle différence
                newDifferenceEpsilon = Math.Abs(matrixX[ctrLineB, 0] - calcul);

                //Si aucun écart n'est inscrit dans le vecteur d'écart on assigne la nouvelle différence
                if (matrixXDifference[ctrLineB, 0] == 0)
                    matrixXDifference[ctrLineB, 0] = newDifferenceEpsilon;
                //vérification de la convergence
                if (newDifferenceEpsilon > matrixXDifference[ctrLineB, 0])
                {
                    //On permet au système d'équation de diverger dans une itération mais pas deux fois de suite
                    if (divergent)
                        throw new Exception("La résolution ne semble pas converger.");
                    else
                        divergent = true;
                }
                else
                {
                    //si sa converge on assigne la nouvelle valeur d'écart
                    matrixXDifference[ctrLineB, 0] = newDifferenceEpsilon;
                    divergent = false;
                }
                //On garde en mémoire le plus grand écart entre les résultats d'inconnus pour pouvoir ensuite
                //le comparer au taux de précision (epsilon)
                if (matrixXDifference[ctrLineB, 0] > largestDifferenceEpsilon)
                    largestDifferenceEpsilon = matrixXDifference[ctrLineB, 0];

                matrixXTemp[ctrLineB, 0] = calcul;
            }

            matrixX = matrixXTemp.DeepCopy();

            //Si la précision n'est pas encore atteinte, on recommence le processus
            if (largestDifferenceEpsilon > epsilon)
                matrixX = FindXUsingJacobi(epsilon, divergent, matrixX, matrixXDifference);

            matrixXResolved = matrixX;

            return matrixX;
        }

        public Matrix FindXUsingGaussSeidel(double epsilon, bool divergent = false, Matrix matrixX = null, Matrix matrixXDifference = null)
        {
            double largestDifferenceEpsilon = 0;
            double newDifferenceEpsilon = 0;
            Matrix subMatrixA;
            Matrix subMatrixX;
            double sumElements;
            double calcul;

            //Instanciation de la matrice de résultat au premier incrément
            if (matrixX == null)
                matrixX = new Matrix(matrixB.SizeY, matrixB.SizeX);
            //Instanciation du vecteur d'écart au premier incrément
            if (matrixXDifference == null)
                matrixXDifference = new Matrix(matrixB.SizeY, matrixB.SizeX);

            for (int ctrLineB = 0; ctrLineB < matrixB.SizeY; ctrLineB++)
            {
                subMatrixA = matrixA.SubMatrix(-1, ctrLineB);
                subMatrixX = matrixX.SubMatrix(ctrLineB, -1);
                sumElements = 0;
                calcul = 0;
                //On fait la somme des sous-matrices en éliminant la ligne et la colonne de l'élément courant
                for (int ctrColA = 0; ctrColA < subMatrixA.SizeX; ctrColA++)
                {
                    sumElements += subMatrixA[ctrLineB, ctrColA] * subMatrixX[ctrColA, 0];
                }
                calcul = (matrixB[ctrLineB, 0] - sumElements) / matrixA[ctrLineB, ctrLineB];

                //On calcule la nouvelle différence
                newDifferenceEpsilon = Math.Abs(matrixX[ctrLineB, 0] - calcul);

                //Si aucun écart n'est inscrit dans le vecteur d'écart on assigne la nouvelle différence
                if (matrixXDifference[ctrLineB, 0] == 0)
                    matrixXDifference[ctrLineB, 0] = newDifferenceEpsilon;
                //vérification de la convergence
                if (newDifferenceEpsilon > matrixXDifference[ctrLineB, 0])
                {
                    //On permet au système d'équation de diverger dans une itération mais pas deux fois de suite
                    if (divergent)
                        throw new Exception("La résolution ne semble pas converger.");
                    else
                        divergent = true;
                }
                else
                {
                    //si sa converge on assigne la nouvelle valeur d'écart
                    matrixXDifference[ctrLineB, 0] = newDifferenceEpsilon;
                    divergent = false;
                }
                //On garde en mémoire le plus grand écart entre les résultats d'inconnus pour pouvoir ensuite
                //le comparer au taux de précision (epsilon)
                if (matrixXDifference[ctrLineB, 0] > largestDifferenceEpsilon)
                    largestDifferenceEpsilon = matrixXDifference[ctrLineB, 0];

                matrixX[ctrLineB, 0] = calcul;
            }

            //Si la précision n'est pas encore atteinte, on recommence le processus
            if (largestDifferenceEpsilon > epsilon)
                matrixX = FindXUsingGaussSeidel(epsilon, divergent, matrixX, matrixXDifference);

            matrixXResolved = matrixX;

            return matrixX;
        }

        /**
         * Affichage du système d'équation dans une string selon les matrices définies 
         **/
        public override string ToString()
        {
            string sortie = String.Empty;
            List<string> valuesMatrixX = new List<string>();

            if (!IsResolved)
                for (int ctrCol = 1; ctrCol <= matrixA.SizeX; ctrCol++)
                    valuesMatrixX.Add("X" + ctrCol);
            else
                for (int ctrCol = 0; ctrCol < matrixA.SizeX; ctrCol++)
                    valuesMatrixX.Add(Math.Round(matrixXResolved[ctrCol, 0], 3).ToString());

            for (int ctrLine = 0; ctrLine < matrixA.SizeY; ctrLine++)
            {
                for (int ctrCol = 0; ctrCol < matrixA.SizeX; ctrCol++)
                {
                    if (ctrCol == matrixA.SizeX - 1)
                        sortie += Math.Round(matrixA[ctrLine, ctrCol], 3) + "*" + valuesMatrixX[ctrCol] + " = " + matrixB[ctrLine, 0] + "\r\n";
                    else
                        sortie += Math.Round(matrixA[ctrLine, ctrCol], 3) + "*" + valuesMatrixX[ctrCol] + "+";
                }
            }

            return sortie;
        }

        #endregion

    }
}
